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Wiskunde en Wetenschappen ondersteunen met ICT 


Numeriek, grafisch, symbolisch 
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KU LEUVEN 


16 de T 3 Europe Symposium 


T‘ VLAANDEREN 

www.t3vlaanderen.be 


19-20 augustus 2013 
KHBO Campus Oostende 


Link: http://www.t3vlaanderen.be/fileadmin/media/symposia/symposium2013/ 

Syllabus symposiuml6.pdf 

Growth Process described and solved by Technology - An Example of a Sustainable 
Learning Spiral (Helmut Heugl) 

Random Walks on Graphs, Matrices Economics and Monkey Business (Hans Schnee- 
beli); Fermi Problems (Hans Schneebeli); Is the Calculus a Must in General Education 
(Hans Schneebeli) 

Without Data there is no Science (Ian Galloway) and others. 


femande@quimica.unlp.edu.ar 

Dear Derivians, 

I am attaching a short dfw file that shows a problem that I have calculating integrals that 
should give gamma functions. Is there any way to overcome this undesirable output? 

Greetings Francisco 

#1: a :e Real (0, 


f 5/2 

#2: J EXP (- a-x ) dx 
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2/5 
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#3: 



4 5/2^, 

-x dx 
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00 
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Liebe bUG-Mitglieder, 

Sie werden sich vielleicht fragen, warum 
die Datei fur diesen bNL so groG gewor- 
den ist? Verantwortlich dafur sind die 
vielen Grafiken im ersten Beitrag. Michel, 
Genevieve und Frederick haben uns ihren 
ACA13-Vortrag als ppt-Prasentation zur 
Verfugung gestellt. Ich habe daraus einen 
Aufsatz gemacht und die vielen Figuren 
aus der Presentation iibernommen. 

Michel und Partner haben erfolgreich eine 
bERIVE-Funktionalitat auf den TI- 
NspireCAS iibertragen und uns dabei inte- 
ressante Einblicke in die Funktionsweise 
eines CAS gegeben. 

Es gibt noch einen zweiten Vortrag von 
ACA13 (Malaga) iiber das Losen von poly- 
nomialen Gleichungssystemen mit TT- 
NspireCAS (Michel Beaudin, Gilles Picard 
und Genevieve Savard). Er wird im nachs- 
ten bNL erscheinen. 

biese beiden Vortrtige brachten mich auf 
eine Idee, die sie auf der letzten Seite 
dieses bNL f inden konnen. 

Erik van Lantschoot hat sich nochmals mit 
der schonen Briicke in seiner Geburtstadt 
beschaftigt. babei zeigt sich sehr deut- 
lich, dass die Geschichte der Mathematik, 
der Ingenieurwissenschaft und die Ge- 
schichte iiberhaupt sehr vernetzt sind. 
Internetquellen mit den digitalisierten 
Texten aus vorigen Jahrhunderten ma- 
chen die Beschtiftigung sehr spannend. 
(Testen Sie Ihre „humanistische Vergan- 
genheit" am lateinischen Text von Cornell's 
Alewyn!) 

Auf der Info-Seite finden Sie einen Link 
zu den Materialien einer erfolgreichen 
Veranstaltung, die von unserem Mitglied 
Guido Herweyers in Oostende, Belgien, 
organisiert wurde. Einige Vortrage und 
Workshops sind in Englisch, daneben gibt 
es Beitrtige in Franzosisch und in Hollan- 
disch. 

Viele GriiGe bis zum nachsten Mai 


bear bUG Members, 

You might wonder about the size of the 
bNL-file. This is due to the many figures 
presented in the first contribution. 
Michel, Frederick and Genevieve provided 
their ACA13 (Malaga) lecture in form of a 
ppt-presentation which I converted to a 
written paper. Many graphs had to be 
copied and pasted. 

Michel and his partners were successful in 
transferring a bERIVE functionality to 
TT-NspireCAS. In their lecture they give 
much insight how a CAS is working inter- 
nally. 

There is a second ACA13 lecture treating 
solving polynomial systems of equations 
with TT-NspireCAS (Michel Beaudin, 
Gilles Picard and Genevieve Savard). It 
will be published din the next bNL. 

Working through both papers I got the 
idea which is presented on the last page 
of this bNL. 

Erik van Lantschoot was busy with the 
historic landmark in his native town once 
more. We can see how history of mathe- 
matics, history of engineering and history 
in general are connected. Internet re- 
sources offering digitized papers from 
earlier centuries made reading the paper 
really exciting. (Test your humanistic edu- 
cation reading Cornell's Alewyn's book 
written in Latin language.) 

The Info page provides a link to materials 
of a very successful conference organized 
by our member Guido Herweyers in Oos- 
tende, Belgium. The keynotes and some 
other lectures are in English, and there 
are lectures in French and butch. 

Best regards until next time 



Download all DTVL-Derive- and TI-f!les from 

http : / / www . austromath . at/ dug/ 
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The DERIVE-NEWSLETTER is the Bulle- 
tin of the DERIVE & CAS- 7/ User Group. 
It is published at least four times a year 
with a content of 40 pages minimum. The 
goals of the DNL are to enable the ex- 
change of experiences made with DERIVE, 
TI- CAS and other CAS as well to create a 
group to discuss the possibilities of new 
methodical and didactical manners in 
teaching mathematics. 


Editor: Mag. Josef Bohm 
D'Lust 1, A-3042 Wiirmla 
Austria 

Phone: ++43-(0)660 3136365 

e-mail: nojo.boehm@pgv.at 


I II 

Contributions: 

Please send all contributions to the Editor. 
Non-English speakers are encouraged to 
write their contributions in English to rein- 
force the international touch of the DNL. It 
must be said, though, that non-English 
articles will be warmly welcomed nonethe- 
less. Your contributions will be edited but 
not assessed. By submitting articles the 
author gives his consent for reprinting it in 
the DNL. The more contributions you will 
send, the more lively and richer in contents 
the DERIVE & CAS- 77 Newsletter will be. 


Next issue: December 2013 


Preview: Contributions waiting to be published 

Some simulations of Random Experiments, J. Bohm, AUT, Lorenz Kopp, GER 
Wonderful World of Pedal Curves, J. Bohm, AUT 
Tools for 3D-Problems, P. Luke-Rosendahl, GER 
Hill-Encryption, J. Bohm, AUT 

Simulating a Graphing Calculator in DERIVE, J. Bohm, AUT 

Do you know this? Cabri & CAS on PC and Handheld, W. Wegscheider, AUT 

An Interesting Problem with a Triangle, Steiner Point, P. Luke-Rosendahl, GER 

Graphics World, Currency Change, P. Charland, CAN 

Cubics, Quartics - Interesting features, T. Koller & J. Bohm, AUT 

Logos of Companies as an Inspiration for Math Teaching 

Exciting Surfaces in the FAZ / Pierre Charland 's Graphics Gallery 

BooleanPlots.mth, P. Schofield, UK 

Old traditional examples for a CAS - what's new? J. Bohm, AUT 

Truth Tables on the TI, M. R. Phillips, USA 

Where oh Where is It? (GPS with CAS), C. & P. Leinbach, USA 

Embroidery Patterns, H. Ludwig, GER 

Mandelbrot and Newton with DERIVE, Roman Hasek, CZK 

Tutorials for the NSpireCAS, G. Herweyers, BEL 

Some Projects with Students, R. Schroder, GER 

Dirac Algebra, Clifford Algebra, D. R. Lunsford, USA 

Treating Differential Equations (M. Beaudin, G. Piccard, Ch. Trottier), CAN 

A New Approach to Taylor Series, D. Oertel, GER 

Cesar Multiplication, G. Schodl, AUT 

Henon & Co; Find your very own Strange Attractor, J. Bohm, AUT 
Rational Hooks, J. Lechner, AUT 

Simulation of Dynamic Systems with various Tools, J. Bohm, AUT 
Space Curves with adjustable Curvature and Torsion, P. Trebisz, GER 
How reliable is mathematical Software? Behrooz Khavari, IRAN 
Two Good Turns a.o., D. Halprin. AUS 
and others 


Impressum: 

Medieninhaber: DERIVE User Group, A-3042 Wiirmla, D'Lust 1, AUSTRIA 
Richtung: Fachzeitschrift 
Herausgeber: Mag. Josef Bohm 
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Integration of Piecewise Continuous 
Functions 


Michel Beaudin, Frederick Henri, 
Genevieve Savard 

ETS, Montreal, Canada 


Abstract 

Piecewise functions are important in applied mathematics and engineering 
students need to deal with them often. In Nspire CAS, templates are an easy 
way to define piecewise functions; in DERIVE , linear combination of indica- 
tor functions can be used. Nspire CAS integrates symbolically any piecewise 
continuous function — and returns, as expected, an everywhere continuous 
antiderivative — as long as this function is not multiplied by another expres- 
sion. DERIVE knows how to integrate sign(a x + b) f{x) where / is an arbi- 
trary function, a and b real numbers and “sign” stands for the signum func- 
tion: this is why products of a piecewise function with any other expression 
can be integrated symbolically. This will be the first part of our talk. 

In the second part of this talk, we will show some implementations that will 
allow Nspire CAS to integrate symbolically products of piecewise functions 
with expressions: the starting point was the discovery of a non-documented 
function of Nspire CAS. Examples of various operations between two pie- 
cewise functions will be given. As a final example, we will show how we 
have defined a Fourier series function in Nspire CAS that performs as well as 
DERIVED built-in “Fourier” function. 


Overview 


ii 


1 . 

1 2 ' 
u. 

4 . 
- 5 . 


Integration of piecewise continuous functions: some problems with Nspire 
No problem with DERIVE ! Why? 

Our solution: Programming new functions in Nspire 
Some applications: Fourier Series and deSolve 
Conclusion 
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Integration of Piecewise Continuous Functions: 
Problems with Nspire 


What Nspire does well : 


/>(*):■ 


-I, x <0 

x, 01 x <2 
2 -x, 2 <x <5 
sin(x),x 25 


Nspire has a nice template that helps the user 
define piecewise continuous functions. 


Done 


fl[x) dr 


-3 


/2(x):=J//(x) dx 


*cos( 7 )+cos( 5 )- — 
2 


Done 



x *0 

0 <x £2 

2 <x <5 


0 

© 

0 

0 



A problem arises when go appears in one of the subdomains: 


* 1 , x <0 

x , 0£x<2 

2-x, 2Zx<5 
sin(.t) p x >5 


dx 


X, 

2 


2 

2 - x- 


x^O 
0<x<>2 

2<x^5 


“1 


-cos (x)+cos ( &)- — ( x> 5 ' 
2 


- 1 , [ -»<x <0 

x, 0 ^x <2 
2-x, 2Zx<S 
sin(x),x >5 


dx 


Nspire can’t compute the 
antiderivative of this function... 


- 1 , -w<x <0 

x r 0 £x <2 
2-x, 2Zx<$ 
sin(x),x ^5 


dx 



- 1 , x <0 




Xj 0 ^x <2 , 

2 -x, 2 ^x <5 

f H 

....nor ol' this function. 

L- 


J 

sm(xj 5 ^x<» 1 


J 


-1, x<0 

x, 0^x<2 
2-x, 2^x<5 
smfx), 5Sx<» 


dx 


3/99 


0 
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A problem occurs when the piecewise function is multiplied by another function 
(even a very simple one). 




/*):■ 


-1, x<0 

x, 0<x<2 
2-x, 2<x<5 
sin(x),x£5 


(/*;(x)- cos(x))dx 


1.24938 


-2 


r 

/• i 

6 \ 


exact 

[fl{x\ cos(x))dx| 


L< 

J-2 / 






COS! 


w- 


-1, x<0 

x, 0<x<2 
2-x, 2<x<5 
sin(x),x£5 


dx 


-2 


X^M' cos(x))d.x 


Nspire can’t find the 
antiderivative. 


COS! 


W- 


-1, x<0 

x, 0^x<2 
2-x, 2^x<5 
sin(x),x£5 


dx 




A problem occurs when we multiply two piecewise defined functions. 



-l,.v <-2 Done 

x 2 , -2<x^l 
- 3 ,.v>l 

\ 


-sin(.r),.v^it Done 

2, x>n 

5 

(/ 7 (.v)-/ 2 (.v))d.v 

J -3 

- 4.85714 

□ 

fi[x)'M {f . M 

1 -sinfvj,.v^7t 1 
U 2 , .t>7tj 

i 

-1, .v<-2 1 

x 2 ,-2<x£1 
- 3 , jt>1 J 



Nspire can’t compute the exact value... 
because it does not simplify* the product into a 
^single piecewise function. 
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No Problem with DERIVE. Why? 


• Defining piecewise functions with some built-in functions (CHI, SIGN, STEP) 

• A very useful integration rule 

• Instead of templates, we may use indicator functions to define piecewise functions in DERIVE. 

• In DERIVE , Indicator (CHI), Signum (SIGN) and Heaviside (STEP) functions are built-in; in 
Nspire CAS, only sign is implemented. 


DERIVE uses the following definitions: 



1, x>0 

“y 

SKNCx) 

SIGN(x) = < 

-1, x<0 



/ 

j+ 

* 

ll 

o 

1 1 * 


STEPW= i ± sip) 

CHI(a, x, b) = STEP(x - a) - STEP(x - b) 



Even though STEP and CHI 
are not built-in in Nspire, 
one can easily define them. 


^ ^ } 1+signLv! 

Define step^x}^ — 


D efine du\ a,\\ b)=step{.\-a} -stepi.x~b) 


l- 

SiGKis} 

X 

: 

I i * 



■ y STEFCvi 


We define piecewise functions in DERIVE as a combination of CHI functions. 

For example, if we need the piece of 7(x) between -2 and 1, we just multiply /(*) by CHI(-2, x, 1) 


For example, if we need the piece of f(x) 
between -2 and 1, we just multiply/^) by 
CHI(-2, x, 1): 


Values at the extremities of the subintervals 
are irrelevant, as far as integration is con- 
cerned. 



-7 -6 -5 -4 -3 -2 -1 1 1 2 3 4 5 
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Another example: 

AO) = 

fl(x) := - x(-“. X, 0) + x-xCO, x, 2) + (2 - x)-x(2, x, 5) + SIN(x).x(5, x, «o) 


-1 

x <0 

-l-x(-oo,x, 0) 

X 

0 < x < 2 

+X-x(0,x, 2) 

2 -x 

2 < x < 5 - 

+ (2 -x)-x(2,x,5) 

,sin(x) 

x > 5 

+sin(x)-x(S, x, co) 


Nspire 


DERIVE 


A(x): = 


-lj x<0 
x, 0<x<2 
2~x, 2<x<5 

sin(x)jX>5 


fl(x) - x(-“. x, 0) + x-x(0, x, 2) + 

(2 - x). X (2, x, 5) + SINCx). X (5, x, «) 


//(x)dx u 

"3 -cos(7)+cos(5) 

2 


11 


- C0S(7) + COS (5) - 


J/*M 


■X, 

„2 


2- x-—. 

2 


x<0 

0<x<2 

2<x<5 


SIGN(x -5). 


COS(x) COS(5) 


x 5 

+ x + — 

4 4 ) 


x * | x - 2| (x + 2)‘|*| COS(x) 


-cos(x)+cos(5) — ,x>5 
2 


Both systems can integrate the piecewise function./? fo). 


-1, x<0 


Nspire 
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We have seen that Nspire CAS is unable to integrate symbolically a product of a piecewise function 
with another expression. 




-1, x<0 '' 

o 

o 

x , 0<x<2 

2-x, 2<v<5 

1 

J ' 

/ i 

sm\xj,x>5 } t 


Can DERIVE find the antiderivative? 


6 

/ fl(x)-COS(x) dx 

-2 

C0S(12) COS(IO) 

+ C0S(5) - 3 • SIN(5) + 2*COS(2) + SIN(2) - 1 

4 4 


and approximated: 


1.249382009 

The definite integral (above) and the antiderivative (below): 

f fl(x)-COS(x) dx 


SIGN(x - S). 


COS(x) 

2 


+ 


2 

SIN(x) 


4 



•SIN(x) + 


C0S(10) 

8 


C0S(5) 

2 


3-SIN(5) 1 

SIGN(x - 2)* (COS(x) + (x - l)-SIN(x) 

2 8 J 


- C0S(2) - SIN(2)) + SIGN(x) 
2 _ 




SIN(x) 


SIN(x) 


DERIVE is able to compute the antiderivative of//(x)cos(x). 
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We have seen that Nspire CAS is 
unable to integrate symbolically a 
product of 2 piecewise defined 
functions. 


Let’s compute the same integra- 
tion with DERIVE. 




# 1 : 

# 2 : 

#3: 

#4 

#5 

#6 


fl(x) := - x(-». x, -2) ♦ x »x(-2 » x, 1) - 3-xCl. x, ») 



-4.857143715 


fl(x).f2(x) 


#7: SIGNCx - n). 


SIGNCx ♦ 2). 


2 2 
Cx l)»SIN(x) x 1 


SIGNCx - 1). 


SIGNCx ♦ 2). 


SIGNCx - 1). 


4 

2 

Cx 

2 

«• 3). SINCx) 

2 . 

2 

X 

3 


4 

2 

2 . 

2 


2 


Cx 

+ 1). SINCx) 

X 

1 


4 

2 

2 . 

2 


2 


Cx 

♦ 3). SINCx) 

X 

3 


- SINCx) - 2 


2 J 


+ SINCx) - 2 


DERIVE unifies the product into a combination of SIGN functions. 


The reasons: 

1 . Piecewise functions are defined as a combination of CHI functions (and this simplifies to 
SIGN functions). 

2. DERIVE knows the rule: 

See DERIVE ’s “stepwise simplification” showing the rule: 

#13: J f(x)*SIGN(a*x + b) dx 


J F(x) «SIGN(a-x + b) dx SIGN(a*x + b) • 

r 

J F(x) dx - SUB ST 

' b y 

J F(x) dx, x, 



k. a 
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This rule is combined with the following rule: 
#15: | x | 

| x I -> x • SIGN(x) 


when DERIVE computes an integral involving an absolute value. For example, 


f I 2 

J x + 5*x + 


dx 


SIGNCx + 3) • 


3 2 

(2*x + 15- x + 36- x + 28). SIGNCx + 2) 



Our Solution: Programming New Functions in Nspire 

Because Nspire CAS is able to integrate symbolically a unique piecewise function (as long as no infin- 
ity appears in the domain!) 
we thought : 

a) to transform the product fl (x)-/2(x) into a single piecewise function, 

b) to “remove” every occurrence of “infinity” in the domain, 

c) and to use the built-in integrator to compute definite or indefinite integrals. 


own - attractive - templates instead of using indi 


We want Nspire CAS to continue using its 
cator functions. 

Our colleague Frederick Henri 
(“Fred”) has programmed some 
simple but quite efficient 
functions. 



a) grouper fct groups in a 
single piecewise function an 
expression that contains one or 
more piecewise subexpressions. 

b) fct_sans_inf ini removes 
every occurence of oo or -oo in 
the domain. 

C ) integral_mcx symbolically 

integrates (indefinite integral) 
and integral_mcx_d 

computes exactly the definite 
integral using the built-in 
integrator. 
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A few Examples: 



We observe the same result when using DERIVE. 


#1: fredl(x) := 2.x(-“, x, 0) + x(0, x, «) 

#2: fred2(x) := - x*x(-«, x, -10) + 2*x«x(-10, x, 5) + SIN(x)*x(5, x, ») 



#3: 


TT 


S 

-l 


fredl(x)*fred2(x) dx 


#4: 




2 

TT - 2 


kit_ets _Jh \groupevJ r ct[ft'edl{x ) • fi'ed2[x),x\ 

-2 x, -~<x<-10 
2- x, 0<x<5 
4- x, -10^x<0 
sin(x), 5<x<°° 


71 

{/red t(x) • Jred2(x))dx 
-1 


7.8696 

kit_ets _Jh \integral_mcx_d[ft'ed l[x) • ft 'ed2[x) t x t - 1,tt) 

71^-2 

( 

71 ^- 2 ) ►Decimal 


7.8696 
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grouper_f ct also works with exponentiation: 



Let’s explain some simple algorithms and show some code. 

First of all, in order to manipulate piecewise functions, we need a command to extract pieces of the 
function. Extraction is not documented into Nspire CAS user guide, but the following “discovery” 
saved us! 

part () 



UNIVERS. 

r\ nspire 


Le site dedie a la communaute T1-Nspire n 


Maths & 
Sciences 



ACCUEIL 

REUSSIR AU lycEe 

REUSSIR EN PREPA 

REUSSIR AU CAPES 

FORUM 

FAQ 

ZONE MEMBRE 









Manuels 


Le site dedie & (a communaute Tl-Nspin 


Index Sujets recents RechercN 


Bienvenue, Invite 
Nom d'utilisateur [" 


Exemples d'utilisatlon 

Premiers Pas 
Exercices 

Programmes Utiles 


Mot de passe perdu ? Nom d'utilisatd Telechargements 


1 

souvenir de moi r | connexion | 


Connectez-Vous 


Nom d'utilisateur 

r 

Mot de passe 


■ Umvers Tl-Nspire le forum ► Le coin aes aeveioppeurs » Discussions generates sur le d6veloppement 

1 □ Acceder aux parametres d'une fonction definie par (1 visualisation(s)) (1) Invite 


E 


REPONDRE AU SUJET 


Page: | 


SUJET: Acceder aux parametres d'une fonction definie par 

■ Acceder aux parametres d'une fonction definie par iiya 4 mois, i semaine #5252 



Bonjour, 

Je cherche une fagon pour une programme en Tl-Basic (Nspire CX CAS) d'acc6der 
aux parametres d’une fonction pre-existante definie par morceaux 

nombre de conditions (ou de morceaux) 

conditionl . valeur de la fonction si conditionl est vraie 

condition2, valeur de la fonction si condition2 est vraie 


Par exemple, si f(x) = x+1 , si 2<x<5 


Se souvenir de moi |“ 

Se connecter 


■f Connexion avec Facebook 


• S'enregistrer sur ce site 

• Mot de passe oublie? 


Tl PlanQt 


La lere Tl-Nspire CM-C CAS 
arrive en France - 
deballage' 
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■ Re: Acceder aux parametres d’une fonction definie par iiy a 3 mois, 4 semaines 


Christian 



Bonjour Genevieve, 

II me semble que I'instruction part peut faire ton bonheur. 


HORS LIGNE 


Moderateur 

eeeeee 

Message: 372 
! - 


/<*):« 

x+l,x>2 
2- x, Else 

Termine 

part^x)) 

4 

panM’f). 1 ) 

x+1 

part(/M.2) 

x>2 


2 x 

p art M x )> 4 ) 

x<2 



Re: Acceder aux parametres d'une fonction definie par i 


semaines 


Christian 



Bonjour, 

Ravi que mon information ait franchi I'ocean pour rechauffer le rude hiver quebecois. 
En fait, c 'eta it deja une instruction non documents dans la voyage 200 et la TI-92. 
Philippe Fortin, que tu connais peut-etre, en avait parle dans une brochure sur le 
calcul formel. J'ai eu I'occasion une ou deux fois de I'utiliser. 

Apres, c'est affaire de memoire: la mienne n'est pas tres bonne, mais part n'est pas 
bien difficile a memoriser. 

Je serai interesse par ton programme, Genevieve. Reviens nous en parler sur le 
forum, quand il sera pret. 

Bien amicalement, depuis la Normandie, a 0°, comme Paris, 

Christian 


h:=Sx+9yz 

8 .v+9 yz 1 

part(/i,o) 

M_|_ II 

part (//J 

2 

pait (/i, 1 ) 

8 - x 

partf//, 2) 

9 yz 


+ 

8x qyz 


partf//, 2 ) 

9 yz 

part (9 -yz) 

2 

partly z,o) 

II . II 

part(9_vz,l) 

9y 

part(9_yz,2) 

z 

+ 


* 

* 

00 

At 

* 

z 

9 

y 
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part ( ) is similar to Maple’s op ( ) 


unifier fcts (fl (x ) , f2 (x) , x) 


fll *) 

fJ2 

f 13 (-1 x < —2 f/11 dom(/ll) 

-<r> 

, ^ 

: -i 

— J ? J Z S x S l aomyj izj 

1-3 x > 1 (/13 dom(/13) 

5 .4 -3 

V J 

: -i 

1 1 1 i If 2(x)=j S1 " w = Km w 

>- 2 x > tt lf22 dom(f 22) 

fll/pl 

-3 -4 -3 

i , , , 

-2 i | i 2 3 4 3 

dom(fll) fl dom{f 21) n domain(J 11/ / 21) 



f f 12 

N 

P3 

-3 -4 -3 

/H 

: -i 


i 2 3 4 3 J 

/22 x 

-3 -4 -3 

: -i 

V 

J 

12 3 4 5 


fl2/J21 1 


•3 -4 -3 

-$ -i 

i 

i 2 y 4 3 


dom(J 11) 0 dom(f 21) n domain(f 11/ / 21) 
dom(f 12) n c/ot?i(/ 21) fl domain{f 12/ / 21) 


/J/ /J2 


^1 

* 

-J -4 -3 -2 4 

f 21 . 


E j J 

4 3 

J? 

-5 -4 -3 -2 -1 1 

L 2 3 

L J 

4 3 


/i3#27 dom(/ll) fl dom(J21) n domain(J 11 // 21) 

J dom(/12) 0 dom(J21 ) n domain(/12//21) 

dom(/13) n dom(J2\) n domain(J\2/f2 \ ) 
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fll 1 { I2 

i P \ 

‘1 

-3 -4 -3 -2 -1 

i 2 3 

4 ! 

/22 x 

■5 -4 -3 -2 -1 

1 2 3 

1 , , 

4 3 

•5 -4 -5 -i 

i 2 

k y 

r 4 

fll/pl /12/J21 

J13/J21 fl3/J22 x 

-3 -4 -3 -2 -1 

i 2 3 

4 3 


/5(.v) : =kit_ets _Jh Wi ifia_fcts [fl[x)" I" J^x), .v) T ermine 


A>J 


sin(.v; 


sin(.v) 

3 




- t .x-nlO ti*0 and v<-2 

,x-nir n*Q and -2<r£l 
t x-nl2 and 1 <jr£it 


sml.v 


— , 7t<V<°° 

'y 


dom(Jll ) fi dom(J 21) fi domain(J 1 1//21) 
rfom(/12) H dom(/21) fi domain(J 12//21) 
dom(/13) n do?n(/ 21) n domain(/ 13// 21) 

dom(/ll)nrfom(f22)=0 
dom(/12) H dom(f22)=0 
dom(f 13) n do?n(/22) n doma[n:(/13//22) 


Then fct sans inf ini ( f 3 (x) ) 


kit_ets Jh 'fct_sans_mfin 


1 


sin(.v) 


,.x-nl3 • 7t*0 and .v<-2 


-.v 


sin(.v) 

3 


— ,x~nl4 • 7t*0 and -2<x&\ 
,x~nl5 • 7t*0 and 1 < v^Tt 


sin(.v) 


— , .V>71 


I 


grouper_f ct (f , x) : generalizes unif ier_f cts by working recursively (in case of 
more complicated functions). 


g(jt):=(6* jc+l}' 

1 f 
5+ 

i 

i. t 

0, jc<3 

1 

Done 

j 

f i- ■ 

grouper _JcAg\x , 

■ ^ 

I x 1 

I I I 1 

< 

5* |6* Jt+l), -°°<jr<3 

r ' l' 9 'l 

\ 6 * x+ 1 j * \x " +5 1 , 3<x< 00 
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Program grouper_f ct: 


grouper^fct (f , x) : = 

operator part(f, 0) 

If f doesn't contain a piecewise subexpression Then 
Return f 

Else 


f 1 (x) : = grouper_f ct (£>art (f , lT|, x) 
f 2 (x) : = grouper_f ct {part ( r , 27 , x> 

Return unif ier_f cts (fl (x) , operator, f2(x), x) 
Endif 


{ 

1 

0, jr<3 V 

5+ 


1 l 

v“..v^3 /, 



grouper_fct (£, x) := 

operator part(f, 0) 

If f doesn't contain a piecewise subexpression Then 
Return f 

Else 

f 1 { x ) : = Igrouper fct(part(f, 1), x) 


f 2 { x } : = grouper_fct (part ( f , 2), x) 

Return unif ier_f cts ( fl (x) , operator, f2 (x) , x) 
Endif 


(6- -t+l) 


H 


0 , x< 5 ' 
x^&Sl 



gr oupe r_f c t { f , x ) : = 

operator := part(f, 0) 

If f doesn't contain a piecewise subexpression Then 
Return f 

Else 

f 1 (x) : = grouper_fct (part (f , 1), x) 
f2 (x) grouper_f ct {part ft, 27], x) 

Return unifier_fcts (f 1 (x) , operator, f2 (x) , x) 
Endif 



(6- .v+l)- |s+| 

0, v<3 \| 

6jc+1 

5+ 1 

0, x<3 
x 2 ,x>3 




fl (x) 

part(f f 2) 
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grouper_fct (f , x) :- 

operator := part(f, 0) 

If f doesn't contain a piecewise subexpression Then 
Return f 

Else 


fl(x):= grouper_fct (part (f , 1), x) 
f2(x):= grouper_f ct (part (f , 2), x) 



Endif 


(.6-v+l)- 


5+ 

0, ,v<3\ 1 




5- (6- .t+l), -°°<v<3 

(6- .v+l) • (v"+5),3Sjr<°° 
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Some Applications: Fourier Series and deSolve with Nspire 

Let us recall that if an expression / of the variable t is defined over the interval t\ < t < tl and ex- 
tended outside the interval by periodicity (the period being P = tl - 1\), then the Fourier polynomial of 
order n of/is the following trigonometric polynomial: 


1 Cj xM 2 

\ f dt+y f cos 

ffAti-ti/ 


Iknt 

tl - tl 


dt cos 


IkTit 


that is: ^ 


+ ^a k cos 

k = 1 


V tl-tly 
r Iknt ' 


+ - 


tl-t\ 


a / 

J/sinl 


+ b k sin 


n 

(”)) 


2 k7Tt ^ 7 . ( 2 k7Tt 

dt sm 


(1 - A 


And this is DERIVED definition from the library “Applications of Integration”: 


1 t2 

FOURIER Cy# t, tl, t2. n) 5 y dt 


tl - tl tl 
2 


tl - tl m_-l 


t? 


t2 


i' 2'Tr-rru-t ‘W ( 2 -it ■ itw_- t ^ ^ f ( 

COS ■ yCOSl dt + SIhI U y-SIH dt 

V tl - tl J l t2 - tl J l t2 - tl J J t t2 - tl J 


tl 


tl 


At ETS, when students need to compute the Fourier coefficients of a periodic signal, they use their 
TI-Nspire CX CAS handheld to compute the integrals (for the Fourier coefficients). 

Then, they store the coefficients and are able to produce any partial sum in exact arithmetic. 


Here is an example. Students are asked to find the Fourier 
series of the following 2 7U - periodic signal. 

„ [ 7T, 0 <X<7T 

./(*)= ' ^ ./(* + 2*) = ./(*) 
-a- + 2;z\ n<x<2n 


The signal is neither odd nor 
even. So, using Nspire CAS, we 
compute the Fourier coefficients, 
splitting the integrals ourselves 
into two parts! 


/>(*):=■ 

it , Q<x<ti Done 

2’ 7 t-x,n<x£ 2 - n 

1 P 

(70:=— 

it 

2 71 3 ' Tl 

fl[x) dx 2 

0 

«(»):=-■ 

71 

7t 

(tt cos(/J r x))dx+ 

Jo 

2' 7! 

((2 7t-x)' COs(tf’ x))dx 

71 

< 

f . 

We split the 

integrals into 

two parts. 

C J 

Done 

j 

*w=- 

51 

7T 

(jr sin(/r -x))d.\ * 

Jo 

”2 * 1 

((2 sin(/j x))dx 

71 J 


Done 

dm) J 

^Then we “inform” Nspire CAS 
that “n” is an integer (in order to 
^simplify the Fourier coefficients), j 

{■if 1 -! 

Ttl^ IT 

b(m) 

1 

nl 
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As mathematics teachers, we are comfortable with this procedure and don’t see any reason to stop 
using it. But on the computer algebra side, being able to automate this procedure is something interest- 
ing. 

We are still teaching some integration techniques despite the fact that the CAS has a built-in integra- 
tor! So, why not define, in Nspire CAS, a “Fourier” function like the one DERIVE has? 

The main goal is to have a Fourier series function able to work in exact mode for piecewise signals. 
This is where the function integral_mcx_d will be useful, replacing the TI’s built-in integrator. 

So, we have defined a “Fourier series function” in Nspire CAS. Using the same syntax as DERIVE and 
replacing the built-in TI integrator by the integral_mcx_d function. 


Let’s check the result: 

#1: f(x) := ji *x(0i x, n) + (2 -tt - x)>x(n, x, 2 -tt) 

#2: FOURIER(f(x), x, 0, 2-n, 5) 

2*COS(5 *x) SIN(5*x) SIN(4-x) 2.COS(3.x) SIN(3.x) SIN(2.x) 2.C0S(x) 3 -tt 

#3: + + + + + SIN(x) + 

25./T 5 4 9-jt 3 2 tt 4 


/*):- 

71, 0<X<71 

2‘ 7I-X,7l<X<2' 71 

Done 

0 \ 


kit_ets_mb\fourie7{fl{x\,x t 0, 2 • tt, 5) 




-2- cos(5- x) ( sin(5- x) ( sin(4- x) 

2- cos(3- x) t sin(3- x) t sin(2- x) 2- cos(x) t ( ^ ( 3- tt 



25-tt 5 4 

+ + +sin^xj+ 

9-tt 3 2 tt 4 



p 20 


B-H-S: Integration of Piecewise Continuous Functions 


D-N-L#91 


Let us mention another interesting application. 

With the command “deSolve”, Nspire CAS can easily solve second order differential equations with 
constant coefficients ... as long as the RHS of the differential equation consists of a single piece. 

Nspire CAS solves a linear second order ODE by using the method of variation of parameters. This 
method involves computing integrals. 


For example, if we try to find a general solution to 


y' + 5y' + 6y = f(t) 


where 



/<0 

/>0 


We will find this: 


sin(r),f<0 Done 


r e',t> 0 

r 




< 


N 


deSolve(y"+5- y+ 6- y=M>*»y) -3. / 

y--e 

H 

sin(r), t<0 \ 
/• e~ ! ,t>0 j 

di 

+e" 2 ' '• 

H 

O O 
V At 

't * 

V) ’C. 

|Jd/- 

c2 e 2 ; +ci e" 3 ; 


V 


) 




) 



Observe the two integrals that Nspire CAS can’t compute. 



When the RHS of an ODE is piecewise, Laplace transforms are usually used. But in this example, 
t can accept negative values. 


With its command “DSOLVE2”, DERIVE can solve the former ODE because DERIVE can compute 
the last integrals. 


So, we have programmed our own desolve2_gen command. 

Our method is still using the variation of parameters but the built-in integrator of Nspire CAS is re- 
placed by Fred’s function integral_mcx_d. 

Let us do the example once more: 


deSolve(y"+5 ■ y+ 6 • y=j(t), ty) 

y=-e~ 3 ' '■ fe 3 'd 

O 

V 

. W 

d,+e' 2 '- f e 2 '•[ 

sin(/), /<0 \ 

dH-c-J e' 2 4c 1- e’ 3 ’ ' 


1 1 

a 

IV 

0 

1 l 

iv 

0 



/(*= 


sin(/),/<0 
/ e' 1 ,^ 0 


Done 




kit_e ts_m b \de solve 2_ge /? ( 5, f) 


This command is able to find a general 
solution to the ODE y” + 5y + 6y 


cl e' 2 *-c2 e ' 3 ‘ 1 C ° S ^ i Sin ^ 
10 10 


g" 3 ' *• ( 5 - (2- /-3)- e 2 ' t+4- (5- ci+4) e ; -20- c2- 3 ) 


20 


»<f<0 


,0</<°° 


rep\- 


-2- 1 -3 / cost/) sin(/J 

ci e -c2 e * 1 — + — — , 


-co</<0 


10 10 

' 3 ‘ *• ( 5 - (2- /-3)- e 2 ' ; +4 ( 5 - ci+4)- e / -20- c2- 3 ) 


0</<~ 


20 


It’s always important to verify the answer! 


-2 / -3- / cosh) sin(/) 

cl e A -c2 e 0 1 — + — — , 

10 10 

e~ 3 ' *■ ( 5 - (2- /-3)- e 2 ' r +4 (5- ci+4) e r -20- c2- 3 ) 


-«><f<0 


20 


# 0</<°° 


kit_ets _Jh\grouper_fct 


[rep)+5- — (rep)+6- rep,t 

\dt 2 dt 


sin(/),-“<?£0 
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Conclusion 

• At ETS, we have adopted TI technology. It started in 1999 with the TI-92 Plus; 
then V200 in 2002 and Nspire CAS CX since in 2011. 

• The CAS system is appropriate for engineering mathematics at the under- 
graduate level. 

• But many mathematics teachers are still using CAS software like Maple or 
DERIVE ; as a consequence, they often want Nspire CAS to be able to perform 
as well as these systems! 

• In the past two years, we have been asking TI to launch a new OS version of 
their CAS system. 

• For the moment, most of their efforts have been on the side of the overall in- 
terface of Nspire CAS. 

• This is correct. But, for mathematicians, the CAS engine should be ranked 
first. 

• “If you want something done right, do it yourself’. 

• The mathematicians needed the help of a programmer. Frederick Henri 
started to work with us. 

• A new team of researchers was formed. In this talk, we showed some results of 
our collaboration. 

• With these new functions programmed by Frederick Henri, the built-in inte- 
grator of Nspire CAS can now be extended to products of piecewise functions. 

• One consequence we showed was the definition of a “Fourier function” similar 
to DERIVE’ s one. 

• And when the TI built-in integrator will be able to integrate symbolically 
piecewise expressions, its “deSolve” command will become better. 


The tns file “KitETSFH” is available for download at 
http://www.seg.etsmtl.ca/nspire/COURS/Kit_ETS_FH.tns 
The Website http://seg-apps.etsmtl.ca/nspire/ 


also contains many examples for using TI-Nspire CAS for undergraduate engineering studies. 
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Brussels Gate in Dendermonde, Belgium, revisited 

Dr.-Ing. E. Lantschoot, Lantschoot_Weisel@web.de 

The topic ’’Brussels Gate” (DNL #89) keeps me busy. Lirst of all, I wanted to know the names of the 
military architects who built these three gates (Brussels, Mecheln and Ghent, the latter completey 
dismantled in 1916). An excellent source was the Internet site which provided the name of captain- 
engineer Comelis Alewyn ( https://inventaris.onroerenderfgoed.be/dibe/relict/48864) . 

Then, I tried to find out more about this person: curriculum, publications, etc. Alewyn lived between 
1788 and 1839; in 1808 he published, in Latin, a treatise on curved lines ( De lineis spiralibus, 
Tentamen academicum ), substantiating his expertise in mathematics c * ] . No other publications were 
found. 

Therefore, I presumed that captain Alewyn did not invent this type of drawbridge, that he relied on 
existing techniques, but where did those techniques come from? 

Having some knowledge about the history of civil and military engineering, I felt that Trance must 
have been the source. Indeed, french military engineering features, in the footsteps of Vauban (1633- 
1707), an uninterrupted succession of engineers (to name a few: Belidor (1697-1761), Monge (1746- 
1818), Poncelet (1788-1867), and others) who were not only prominent in their achievements, but who 
developed new designs and applied the most recent methods of calculus to them. In addition to that, 
Holland had been under the rule of Napoleon during the period 1795-1813, and it can be assumed that 
dutch professionals like captain-engineer Alewyn were familiar with the Trench publications on the 
subject of their activities. 

A drawbridge in Trench is a ”pont-levis”. The Trench Wikipedia site http://fr.wikipedia.org/wiki/Pont- 
levis provides an overview of all different types of drawbridges, and classifies them according to their 
system of counterweight. The design of the Brussels gate is found in § 2.6, ” constant counterweight on 
a guiding rail ” and it is attributed to Bernard de Belidor (1697-1761). Note: the Wikipedia site must 
be read in Trench; had we selected another language, we would have obtained one more discussion of 
the topic in the new language, not a translation of the Trench original. The publications of Belidor can 
be consulted in digitalized form in http://gallica.bnf.fr/ and his derivation of the drawbridge curve can 
be found in ” La Science des Ingenieurs ”, Paris 1729, book 4, chapt. 5. He called it a Sinusoid, 
although by the nomenclature of today it is an Oval of Descartes, or a subset thereof: Snail of Pascal 
or Cardioid. Another source for the derivation of the Belidor curve is a publication of J.-V. Poncelet, 
’’Traite de la construction des ponts-levis”, Paris 1845, which can be consulted in digitalized form 
under http ://www. e-rara. ch/doi/ 10.3931 /e-rara-3 6 ; the Belidor curve is discussed on pdf-page 18, the 
figure appears on p. 64. An interesting sideline: where Belidor in 1729 works with ’’livres” (pounds) 
and ’’pieds” (feet), Poncelet complies with the (in the meantime, 1795, legalized) metric system. By 
the way, Poncelet is the man who invented the projective geometry, when in captivity in Russia during 
the campaign of Napoleon, 

The Belidor counterweight system was applied occasionally (but without much success) to street-wide 
bridges in the USA; examples can be seen in http://bridgehunter.com/category/tag/belidor-bascule/ or 
in http://historicbridges.org/ , search for ’’Belidor”. 

[1] You can download the Latin book from http ://www. e-rara. ch/zut/content/titleinfo/ 1028122 . (Chapters are: 
Tangentis Applicationes , De Circulo Osculatorio, De Spiralium Quadratura, ...)I found another book written by 
Alewyn: Beschrijving van een ontwerp van sluizen met gekoppelde deuren, welke bij alle waterstanden geheel of 
gedeeltelijk geopend en wederom gesloten kunnen worden. Josef 
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The Descartes Oval 

This curve describes the movement of the counterweight in a Belidor drawbridge system. We are 
using the notations of R. Ferreol in http://mathcurve.com/courbes2d/pontlevis/pontlevis. shtml . [ * ] 

^ p : weight of the bridge [N] 

t q : counterweight [N] 

t a : length of the bridge [m] 

^ b : distance between the plane of the bridge and the center of gravity [m] 

t h : distance of the street level from the deflection sheave O which is the origin of the 
coordinates [m] 

^ 1 : total length of the rope [m] 

> fr : free rope = lenght of rope protruding at the city side when bridge is down [m]. Is also the 

distance OS. 

^ rho, 0 : polar coordinates of the counterweight = point B on the Belidor curve [m, rad]. In contrast 
with R. Ferreol we define 0 in the usual way, i. e. from the x-axis (0 = 0) and turning to the 
left. 

^ a : angle OHT between the bridge and the gate [rad] 

^ O : Deflection sheave 

> H : hinge of the bridge 

► T : tip of the bridge 

^ S : highest point of the Belidor curve 

► W : a weight simulating the force of the winch. 

The fundamental consideration is that the equilibrium of the point B on the Belidor curve implies that 
the increases and decreases in potential energy, which arise when rho becomes larger, do have a 
constant sum. 

The graph depicts a most general situation: a and h are not equal, there is some free rope fr at the city 
side when the bridge is down, and a constant additional force w is generated in the rope by the pull of 
a weight w (we will later show to what use). In our mind, we now stretch rho from rho = fr to some 
arbitrary higher value, and observe four quantities: 

1) Point T moves upwards, its potential energy increases by p/2*a*cos(a), 

2) Point B originally was at S and possessed a potential energy of -q*fr, 

3) When B comes down from S to B, there is a decrease in potential energy of q*rho*sin(0). 

4) The length of the second rope from O to W is immaterial; if W originally was at O with 0 potential 
energy, it now is at -(rho - fr) *w. 


1 http ://mathcurve .com/ is a great website. I recommend visiting this site (if you are interested in 
curves and surfaces), . . . Josef 
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By stating that the sum of those four contributions equals 0 we obtain the equation 

(1) sumpot = p/2*a*cos(a) - q* (l-V(a A 2+h A 2)) + q*rho*sin(0) - (rho - 1 +V(a A 2+h A 2)) *w = 0. 

This sum is computed, in algebraic notation, by the program oval (). If numerical values are desired, 
we call data() and answer the questions. belidor() finally applies statistical methods to the results. But 
yet back to sumpot. 

Applying the cosine rule to the triangle OHT, we can eliminate cos(a). sumpot now appears in terms 
of a, h, p, q, and 1, which are the design parameters of the engineer. 

sumpot = 0 is a quadratic in rho with 0 as a parameter. However, it is still an equation with dimension 
energy; what we are looking for is an equation with dimension length. We can multiply both sides of 
sumpot = 0 by any non-zero quantity, and -h/(l*p) is a good candidate. We thus obtain 

(2) sumdis = sumpot*h/(l*p) = 0. 

This is a quadratic in rho, with 0 as a parameter. It is a physical implementation of a class of mathe- 
matical curves known as ’’Oval of Descartes”. These curves are very important in mechanical and 
optical engineering, and no one less than the young James Clerk Maxwell started his career with an 
article on them (1846). 

Suppose now, that we list the coefficients of the powers of rho of (2). This is easily done by the 
command polycoeffs(sumdis,rho). This command delivers a three-item list desclst and the Oval of 
Descartes appears as: 

( 3 ) desclst[l] * rho A 2 + desclst[2] * rho A l + desclst[3] * rho A 0 = 0. 

Its two roots rtl and rt2 are found by the command zeros(sumdis, rho). When plotted, the oval of 
Descartes appears as two nested ovals, each one corresponding to one of the roots. The root which 
generates the outer oval is the one which pertains to our problem. 

If we have supplied new data with the program data(), and if we turn to the graph, we can play around 
with the counterweight q, and we will soon observe, that the selected value fr is implemented as 
desired for values of q above the minimum value qmin, but not implemented as desired for values of 
q below qmin. 

Thinking of it as an engineer, it makes no sense to use a counterweight lighter than qmin, because the 
equilibrium could not be achieved. 

But how is this stated in mathematical terms? Suppose rtl is the outer oval. If not, use rt2. We 
suggest that you take the following steps: 

1) rtlJ on the scratchpad. Its formula appears as a numerical expression containing sin(0). 

2) DelVar q <J: The numerical value of q is taken away, and replaced by q as a variable. 

3 ) Now we apply a condition: rtl|0 = -nil and rtl appears in the following form: 

(4) rtl|0=-7r/2 = factor * (|q - qmin| - q + constant). 

If q>qmin, we may forget the || signs, and see that q disappears from the expression. Otherwise, 
q remains in the formula. 
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If we look at the outer oval rtl as a whole (modify the scale of the graph to obtain the whole outer 
oval) it comes to our mind, that the curve has about the shape of a degenerate form of Descartes, viz. 
the Snail of Pascal, which occurs when 


desclst[3] = 0. The Snail of Pascal is usually written in the standardized form 


(5) rhop = 2*1* (1 + 2* (h/1) * (q/p) *sin(0)). 
and the coefficient of sin(0) is usually referred to as e 


(6) e = 2 * (h/1) * (q/p). 

With foresight, the value of e was already calculated by data() and rhop is contained in the list of 
polar curves displayed on the graph, but must be activated by a check sign. We see that both curves fit 
very well in general, but unfortunately rather poorly in the area around S. 

For the particular example of Dendermonde with a = h = 3.4 m, p = 12000 N, fr = 0.9 m and 
q = 8500 N the street is reached at a point G = (4.80, -3.4) m, which is well in accordance with the 
observations on the spot. But, be careful: HG must = 1. 

Any further study of the curves of Descartes and Pascal, although very interesting, would go far 
beyond the scope of this notice. Fortunately many good textbooks on plane curves are available in 
bookshops and libraries. 

Belidor is satisfied 


In the meantime, it appears that the most important design parameter is q, the weight of the 
counterweight. The engineer draws a circle with length 1 around O, sees where it intersects the Belidor 
curve, and then adjusts q so that the intersection comes down to the street level. 

What we will do now is to prepare the Belidor curve for regression analysis . The polar coordinates for 
8 points are converted to two lists rx and ry of cartesian coordinates. To these lists we apply the 
quartic regression and obtain the exact 4th-degree equation in x (*) of the Belidor curve. Then we 
apply the quadratic regression and obtain eq2, the best possible parabolic approach to Belidor.The 
differences are hard to notice, but they exist: in the first half Belidor is under the parabola, in the 
second half above it and then under it again. The check is easy: it suffices to plot 100 * abs(eq4-eq2) 
and the differences can be seen in centimeters. They are, for V(l A 2-h A 2) > x > 0 at worst about 10 cm. 

I am e-mailing a hint to the Belidor drawbridges at Dendermonde to "mathourist.canalblog.com". If 
you want to know what this Internet site has to offer, do not forget to enable your browser to access 
French sites. 
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Alewyn 's problem solved mathematically 

Beautiful as it is mathematically speaking, the Belidor curve is not the proper practical solution. Any 
climatic disturbance, rain, snow, wind would move and shatter the bridge uncontrollably. Therefore, 
Alewyn had to introduce a disequilibrium, which would (a) forcefully push the bridge on the 
bridgehead when in the horizontal, and which (b) would have to keep the bridge forcefully upright 
when in the vertical position. 

To achieve these two goals Alewyn took advantage of 
the fact that (in practice) the center of gravity does not ly 
on the surface of the bridge but rather b m below it, and 
reasoned that the equilibrium of energies would have to 
be obtained when the energy exerted by a winch was 
included in sumpot. If the force w of the winch were not 
present, the brigde would ly flat forcefully, because the 
counterweight would contribute less force than required 
for equilibrium, and if the bridge stood upright, the force 
of the winch would keep it in equilibrium. 

We can simulate the force w of the winch by a weight W 
moving up and down and hanging on a second rope that 
begins at the tip of the bridge, then goes over the 
deflection sheave and moves straight down to W. 

Mathematically speaking, this second rope can be 
treated as an additional force w in the first rope. Assuming that a = h as in Dendermonde, we now call 
data(). The program asks for b (center of gravity below the plane of the bridge) and proposes a 
sensible value for w, viz. the value needed for keeping the bridge in equilibrium when it is in the 
upright position. Higher values are tolerable (to the detriment of the men who power the winches!) but 
lower values would not insure the equilibrium anymore. Looking at the formulas, one sees that w 
plays a role similar to q, the counterweight (**). 



Finally, 

the design invented by Belidor was not very successful. To my esteem, the freedom of the designing 
engineer was too small: there is only one parameter - 1 mean q - he can play with, and even then only 
within rather narrow limits. 


Notes 

(*) It is not immedialtely obvious that Descartes (when stated in Cartesian coordinates) is a 4-th order 
equation in x, but if we plot that equation, we can see that the x-axis is crossed never more than 4 
times. 

(**) The numerical results are approximations only, because we act as if p/2 was still located at the tip 
T of the plane of the bridge, but all we wanted to do was to show the principle of Alewyn’s solution. 

More Links (German): 
http://de.wikipedia.org/Bernard de Belidor 

http://www.deutsches-museum.de/bibliothek/unsere-schaetze/technikgeschichte/belidor/ 



D-N-L#91 


E. van Lantschoot: Brussels Gate in Dendermonde rev. 


P 27 


Define data()= 

Prgm 

:© Supplies or modifies the data 
:DelVar a,b,h,fr,p,q,qmin,w 
:setMode({ 1,4,2,1 }) 

:© Typical values are entered 

:3->a: 4->-h: 0->-fr: 10000 ->p: 0-»b: 0-*w 

:Disp "COMPUTER OR HANDHELD: SET ANGLE UNITS TO RADIANS!" 
:Disp "SET NUMERICAL FORMAT TO FLOAT 4" 

: Request "lenght of bridge = ",a 

: Request "height of deflection sheave above street = ",h 

: 10000 -> p © weight of bridge 

: Request "weight of bridge = ",p 

:If a^h Then: Goto next: Endlf 

:Request "center of gravity below bridge plane = ",b 

:((p*b)/(h))->-w 

:Request "force of winch = ",w 
:Lbl next 

:© provisory value of q = value for equilibrium at 0=-((n)/(2)) 

:q:=((V(a A (2)+h A (2))*p)/(2*h))-w 

:qmin:=q 

:Disp ’’minimal counterweight = ”,q 

:Request ’’length of free rope between 0 and ”&string(h)&” ?”,fr,0 
:If fr<0 or fr>h Then: Disp ’’choice not appropriate”: Stop 
:Else 

:Disp ’’selected length of free rope = ”,fr 
: Endlf 

:l:=V(a A (2)+h A (2))+fr 
:Disp ’’total length of rope = ”,1 
:e:=approx(2*((h)/(l))*((q)/(p))): Disp ”e =”,e 
:EndPrgm 


datd) 

COMPUTER OR HANDHELD: SET ANGLE UNITS TO RADIANS! 
SET NUMERICAL FORMAT TO FLOAT 4 
lenght of bridge = 3 

height of deflection sheave above street = 4 
weight of bridge = 10000 
minimal counterweight = 6250 
selected length of free rope = 1.1 
total length of rope = 6.1 
e =0.82 


Fertig 


i 

1/99 
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Define belidor()= 

Prgm 

:setMode(2,l) 

:rx:=newList(8) 

:ry:=newList(8) 

:0step:=((cos- 1 (((h)/(l))))/(7)) 

:For j,l,8 

: 0 : =((-:i)/(2))+j *0step 

:© select the outer oval, rtl or rt2 

:rx[j]:=rtl*cos(0) 

:ry[]']:=rtl*sin(0) 

:EndFor 
:QuartReg rx,ry 
:eq4 : =stat . RegEqn(x) 

:Disp "exact Belidor eqn = ",eq4 
:QuadReg rx,ry 
: eq2 : =stat . RegEqn(x) 

:Disp "best parabolic approximation = ",eq2 
:EndPrgm A 



Define oval()= 

Prgm 

:set*mode(2,l) 

:Local ant,coa 
:Unlock q 

:DelVar 0,a,b,tr,h,l,p,q,w 

:© cosine rule for OHT. coa stands for cos(a) 

:coa:=zeros((l-rho) A (2)-a A (2)-h A (2)+2*a*h*coa,coa)[l] 
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:sumpot:=((p)/(2))*a*coa+q*(l-V(a A (2)+h A (2)))+q*rho*sin(0)+(l-V(a A (2)+h A (2))-rho)*w 
:Disp "If no interest for Alewyn, disregard the existence of w" 

:Disp "sumpot = " 

:Disp sumpot 

:sumdis:=((“Sumpot*4*h)/(p)) 

:Disp "sumdis =" 

:Disp sumdis 

:desclst:=polyCoeffs(sumdis,rho) 

:Disp "Coefficients of Descartes' equation with 
rho as the variable and 0 as a parameter 
describing both ovals: = desclst" 

:Disp desclst 

:Disp "Inner and outer oval are =" 

:rtl :=zeros(polyEval(desclst,rho),rho) [ 1 ] 

:rt2:=zeros(polyEval(desclst,rho),rho)[2] 

:Disp "rtl = ",rtl: Disp "rt2 = ",rt2 

:© computing an approximation by the snail of Pascal 

:rhop:=2*l*(l+2*((h)/(l))*((q)/(p))*sin(0)) 

:EndPrgm 
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Quasi-Monte-Carlo-Methods 

Josef Bohm, Wiirmla, Austria 


Instead of an abstract 

A quasirandom or low discrepancy sequence, such as the Faure, Halton, Hammersley, Nied- 
erreiter or Sobol sequences, is ’’less random” than a pseudorandom number sequence, but 
more useful for such tasks as approximation of integrals in higher dimensions, and in global 
optimization. This is because low discrepancy sequences tend to sample space ’’more uni- 
formly” than random numbers. Algorithms that use such sequences may have superior con- 
vergence. Faure sequences, in particular, seem to have become popular in mathematical fi- 
nance simulations. [8] 

Introduction 

Monte-Carlo-methods (Monte-Carlo-simulations) are widely used for simulations based on inspecting 
and evaluating many randomly produced samples. The best known example - even in secondary 
school - is approximating n by counting “raindrops” fallen on a sqare with an inscribed quarter of a 
circle. 

MC-methods are used in statistics, physics, medicine, ..., and in financial mathematics (option pricing, 
estimating risks, insurance, ,„), too. You can find a good overview and many links in 

http://en.wikipedia.org/wiki/Monte Carlo method. 

MC-methods use random numbers produced by random number generators which are in fact 
producing pseudo random numbers applying a less or more good algorithm. Many programs need a 
random seed to initiate the agorithm. Same seeds produce identical random number sequences. Some 
programs - like DERIVE - can use the internal CPU-time for starting the process, which gives new 
sequences whenever called. 

Quasi-Monte-Carlo-methods use deterministic sequences instead of random ones. They outperform 
the MC-methods in convergence and error bounds. Both attributes proved to be important in financial 
mathematics applications. 


Once more the traditional MC-method 


I will start reminding on the MC -integration but not 
referring to the - in the meanwhile maybe boring - 
example of estimating n. 

I was inspired by an article in “Spektrum der 
Wissenschaft” - which is the German edition of 
Scientific American - from February 2012 to work 
with QMC-methods. 

The first figure in the journal shows a leaf and how to 
estimate its area. There are 1024 MC-randomly 
distributed points in the square and 41.89% of all 
points are within the area of the l ea f. [Seeer,doftHn,nbutlo,,] 
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I had liked to use the picture of the leaf as background picture in a DERIVE plot window but I was not 
able to reproduce the shape of the leaf by functions. So I decided to take a star for demonstration 
purposes. I defined six inequalities for describing an irregular hexagon, calculated the intersection 
points of the lines and then calculated the area of the polygon. 


(I skip presenting the calculation of the vertices.) 


[il(x, y) := y - 5*(1 - x) > 0 , i2(x, y) := 2«y - 10 + x < 0] 

[i3(x p y) := x - y - 1 < 0, i4(x, y) := 2«y + x - 4 > 0] 

[i5(x p y) := y - 5.(4 - x) < 0 , i6(x p y) := y - 2-x - 2 < 0] 

p(x p y) := (il(x p y) a i2(x p y) a i3(x p y)) v (i4(x p y) a i5(x p y) a i6(x p y)) 


1 


a_tr(P , Q p R) := — 
2 


P - (Q - R ) + Q • (R - P ) + R .(P -Q) 
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Calculation of the area would have been easier if I had used 
TI-N spire (see left below). The result is the same, fortunately! 




TI-Nspire 


GeoGebra 
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To go for sure I plotted the hexagon once more using GeoGebra - the area of the hexagon corresponds 
with the DERIVE - calculated value again. 

I don’t want to loose time complaining about accurracy but proceed by estimating the area of the 
hexagon applying the traditional Monte-Carlo-method using DERIVED random numbers. 

3000 points should be sufficiently enough! 


RANDOM(O) = 2208522468 

mc_polpts := VECT0R( [4- RANDOM (1), 8 . RANDOM(l)] , k, 3000) 


DIM(SELECT(p(v , v ), v p mc_polpts)) = 1412 
1 2 


SOLVE 


32 


3000 

■ a 

1412 y 


= (a = 15.06133333) 


DIM(SELECT(p(v , v ), v, mc_pol pts)) • 32 
1 2 

3000 

DIM(SELECT(p(v p v ), v, mc_pol pts)) . 32 
1 2 

3000 

These are the results of three simulations 
15.2643 - a 


= 15.36 


= 15.22133333 


VECTOR 


• 100 p a p [15.0613, 15.36, 15.2213] 


15.2643 

[1.329900486, 0.6269530866, 0 . 2817030587] 

The maximum error is 1.33%. (Please compare with page 40.) 



Demands on random numbers 

There are three demands on random numbers: 

( 1 ) They must not be predictable. 

(2) They must not be correlated, i.e. they must be independent of another. 

(3) Their distribution must be well-balanced. 

True random numbers fullfil all demands, pseudo random numbers follow a strict - more or less 
complicated - rule, quasi random numbers fullfil only one of the three demands: they are evenly 
distributed. A measure for the balance of the distribution is the so called discrepancy. The discrepancy 
is important for the possible error. You can find much about this on the respective websites. 

The minimised error together with a very much improved convergence made quasi random numbers 
and the corrresponding Quasi-Monte-Carlo-methods a modern tool for financial mathematics. 
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As this is not very elegant, I did once more without using the string operations but conversing, 
reversing and calculating digit for digit: 


corput_(x, b := 2, p, o) := 
Prog 

[o := 0, p := 1/b] 

Loop 

If x = 0 
RETURN o 

o := o + M0D(x, b)«p 
x := (x - M0D(x, b))/b 
P := p/b 


corput_seq(n, b := 2) := VECTOR(corput_(k , b), k, n) 


corput_seq(10) 


■ 1 

1 

3 

1 

5 

3 

7 

1 

9 

5 ' 

1 

. 2 

» 

4 

i 

4 

i 

8 

p 

8 

i 

8 

8 

16 

16 

16 . 


For comparing with DERIVE ' s (pseudo) random numbers I generate 1200 elements of the Corput- 
sequence in order to simulate throwing a die 1200 times and observing the outcomes: 

#8: dice := VECTOR (RANDOM (6) + 1, k p 1200) 

Simplify the expression above, remove the huge expression and proceed counting: 

#9: VECTOR(DIM(SELECT(k = i p k p dice)), i, l p 6) 

#10: [189, 197, 20 5, 212, 198, 199] 

#11: dice_qmc := corput_seq(1200) 


Simplify again - then you can remove the huge expression, proceed counting 




r 

r 

f i 

i + 1 

*v 

\ "N 

#12: 

VECTOR 

DIM 

SELECT 

— < k < 


k, dice_qmc 

,i,0,5 





l 6 

6 

J. 

/ J 


#13: [202, 200, 198, 202, 200, 198] 

Simplifying #8 and #1 1 is important for fixing the sample. Otherwise di ce and di ce_qmc in 
expressions #9 and #12 would be generated new at each call within the VECTOR-command. 

From the Sequence to Quasi Random Points (Hammer sley) 

There are several ways to create sets of random points in the plane, in 3D-space and - and this is the 
interesting part for applications in financial mathematics in many dimensions. 

I start presenting the Hammersley Set described on page 33. 

hammersley(n, lb := [0, 0], rt := [1, 1], b := 2 , a_, b_, x_, y_) ;= 

Prog 

[a_ := rtil - lbil, b_ := rtj,2 - lb|2] 

x_ := [VECTOR(lb|l + i/n-a_, i, 0, n - 1)] 

y_ := [VECT0R(lbi2 + k«b_, k, corput_seq(n , b))] 

APPEND(x_, y_Y 

(lb = left bottom vertex and rt = right top vertex of the rectangle describing the region in the plane. It 
is the uit square by default) 
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[hammersley(lO) , 
0 0.5 ' 

0.1 0.25 

0.2 0.75 

0.3 0.125 

0.4 0.625 

0.5 0.375 

0.6 0.875 

0.7 0.0625 

0.8 0.5625 

L 0.9 0.3125 . 


hammersley(10 , [2, 
2 1.5 

2.1 1.25 

2.2 1.75 

2.3 1.125 

2.4 1.625 

2.5 1.375 

2.6 1.875 

2.7 1.0625 

2.8 1.5625 

2.9 1.3125 


1], [3, 2]), hammersley(10 , [-2, -1], [0, 3])] 
-2 1 


- 1.8 

- 1.6 


-1.4 -0.5 

-1.2 1.5 


-1 


0.5 


-2 - 1.5 


- 0.5 


-0.8 2.5 

-0.6 -0.75 

-0.4 1.25 

-0.2 0.25 

■ 2.5 

2 

1.5 
1 

■ 

■ 

■ 

■ 

0.5 

- 0.5 


1.5 


2.5 


10 quasi random points distributed in three different regions in the plane. (Unit square and base 2 by 
default). Let’s plot 1000 and 2000 points (base 5) now: 


hammersley(lOOO) 



hammersley(2000 , [0, 0], [2, 5], 5) 
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Generalisation: Halton Sequence and Halton Set 


Why only taking base 2? The Halton sequence is an ^-dimensional generalization of the van der Cor- 
put sequence, but instead of using binary representations, a different base {prime number) is used for 
each coordinate. Let’s take p\ = 2 and pi = 2 for x- andy-axis. 

halton_seq(n, b) ;= corput_seq(n , b) 


halton_seq(10 , 2) 


halton_seq(10 , 3) 


hal ton_seq(10 , 5) 


I" 1 1 3 1 5 3 7 1 9 Si 


1 

. 2 

4 

i 

4 

i 

8 

» 

8 

p 

8 

8 

p 

16 

16 

16 

- 

" 1 

2 

1 

4 

7 

2 

5 


8 

1 

10 ' 


i 

. 3 

t 

3 

i 

9 

i 

9 

p 

9 

i 

9 

9 

p 

9 

1 

27 

1 

27 . 


' 1 

2 

3 

4 

1 

6 


11 


16 

21 

2 " 

i 

. S 

i 

5 

5 

i 

5 

i 

25 

25 

i 

25 

p 

25 

25 

25 . 


halton(n, lb, rt, bl, b2, a_, b_, x_, y_) := 

Prog 

[a_ := rt|l - lb|l, b_ := rt|2 - lb|2] 
x_ := [VECTOR(lb|l + k*a_, k, hal ton_seq(n , bl))] 
y_ := [VECT0R(lb|2 + k«b_, k, hal ton_seq(n , b2))] 
APPEND(x_, y_)' 



In http : //planning . c s .uiuc . edu/node2 1 0 .html I found a nice figure comparing the Voronoi diagrams of 
196 Hammersley points and 196 Halton points and of 196 (pseudo) random points. I will try to 
produce the diagrams with 100 points only (Hammersley left and Halton right). 
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I am not so happy with the Halton-Voronoi diagram. 

Then I plotted the Voronoi diagram for 100 (pseudo) 
random points (right) . You can observe the difference 
concerning the “uniform” distribution. 

You may remember Enric Puig’s request on Voronoi 
diagrams in DNL#89. 1 can not give an advice how to 
contruct such a diagram but I can give a hint how to plot it. 

I copied the points to GeoGebra which provides a 
Voronoi[<pointlist>] command and plotted the diagrams as 
you can see above and right. 



Quasi Monte Carlo vs Monte Carlo 

I come back to my introduction talking about approximating n. 


qmcpts := q_m_hal ton([0 , 0], [1, 1], 2, 3, 1000) 


2 2 
x + y = 1 

I 0 0 1 
1 0 
1 1 
0 1 
. 00 . 

2 2 

0<x<lA0<y<lAX + y <1 
qmcpts := halton(1000, [0, 0], [1, 1], 2 , 3) 



r 

'2 2 Y 

DIM 

SELECT 

v + v < 1, v, qmcpts 

l 1 2 )) 


4 -787 

787 

1000 

250 

4*787 

= 3.148 


1000 


I repeated the estimation using 5000 points and 
received 3.1472 (QMC) vs 3.164 (MC). 
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In my next example I calculate the area between;; = 0 and/ix) = sin(2x) +1 for 0 < x < 7i/2. 
f (x) ;= SIN (2-x) + 1 



r 

7T 


sinpts := halton 

2000 p [0 P 0], 

, 2 

. 2 

, 2, 3 

y 


DIM(SELECT(f (v ) > v , v p sinpts)) = 1639 
1 2 

1639- 7T 

= 2.574535179 

2000 

tt/2 

J f(x) dx = 2.570796326 
0 

| 2. 576105975 - 2. 570796326 | 

= 0.002065371319 

2.570796326 



mc_sinpts := VECTOR 


RANDOM (1) , 2 • RANDOM (1) 

L 2 


DIM(SELECT(f (v ) > v , v, mc_sinpts)) = 1634 
1 2 


1634 - n 

= 2.566681197 

2000 

| 2. 566681197 - 2. 570796326 | 

= 0.001600721519 

2.570796326 



Quasi Monte Carlo in 3D Space 

The applications in financial mathematics are very often evaluations of multiple integrals. So we have 
to leave the two dimensions and generalize for more dimensions. We will use the Halton sequences to 
calculate the volume under a surface. 


There are QMC-methods which are better suitable if many dimensions are needed. I will present them 
later, hal ton3d produces a list of n points The region is defined by three vertices of a cuboid 
(lfb = left front bottom, rfb = right front bottom, lbt = left back top). 

halton3d(n p lfb p rf b p lbt p bl p b2, b3 p a_ p b_ p c_ p x_ p y_, z_) ;= 

Prog 

[a_ := rfb|l - lfb|l p b_ := lbt|2 - lfb|2 p c_ := lbt|3 - lfb|3] 

x_ := [VECTOR(lfbil + k«a_ p k p hal ton_seq(n p bl))] 

y_ := [VECT0R(lfb|2 + k«b_, k p hal ton_seq(n p b2))] 

z_ := [VECT0R(lfbi3 + k-c_, k p hal ton_seq(n p b3))] 

APPEND(x_ p y_ p z_)' 


pts3d := halton3d(1000 , [0 P 0, 0], [l p 0, 0] p [0, l p 1], 2 P 3 P 5) 

plot3dl i st(l i st) ;= VECT0R([v] p v p list) 

plot3dlist(pts3d) 



We will estimate the volume below a surface which needs a double integral. Let’s assume that we 
cannot integrate the function. The problem is: 


2 2 

x y 

g(x, y) := l 

2 2 




ptsBd 


hal ton3d(1000 , [0, 0, 0], [1, 


0, 0], [0, 1, 1], 2 , 3, 5) 


DIM(SELECT(g(v , v ) > v p v, pts3d)) = 667 
12 3 


667 


= 0.667 

i nnn 



We count the points whose 3 rd coordinate is 


less than the function value. 


Then we plot the surface and the points under the 
surface (blue) and above (red). 
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The Hexagon from page 32 

At the end of the Corput-Halton-section of my contribution I’d like to come back to the hexagon from 
page 32 and estimate its area with two different Halton sets (different bases): 


polpts := haltort(3GQQ p [0, 0], [4, S] , 2 , 3) 

DIMfSELECTtpO f v ) t v t polpts)) = 1434 
1 2 


32 


SOLVE 


3000 

■ a 

1434 


(a = 15.296) 


15.2643 - 15.296 


15.2643 


0.002076741154 


polpts := halton(3QQQ, 


DIM(SELECr(p(v p v ), 
1 2 


SOLVE 


32 


3000 


a 1433 
15.2643 - 15.2853 
15.2643 


[0, 0] p [4, 8] p 5 , 11) 
v, polpts)) = 1433 


(a = 15.28533333) 


= 0.001375759124 


Errors are 0.2% and 0.1%. Convincing or not? 


Faure Sequences 

When I started writing this contribution and programming the first function I only wanted to refer to 
the article in the Scientific America. Then I was ready presenting the Hammersley set. Curious as I am 
I did some Internet research - 1 shouldn’t have done it. To quote our friend David Halprin: 

“I opened a can of worms.” (in DNL#29 at the occasion of his Super Duper Osculants). I found much 
materials, too many papers - and most of them filled with not so easy mathematics. They wrote about 
van der Corput, Halton, then much about Faure- and Sobol ’ sequences. I came across a famous 
Austrian mathematician who was also busy in this field of mathematics, Edmund Hlawka (Koksma- 
Hlawka Inequality). Finally I was filled with ambition to include the Faure Sequences. They are 
mentioned in all papers treating QMC-methods. 

I must admit that I did not understand the algorithm described in the first paper. I inspected the next 
one, and the next one, ... - surprisingly in most papers I found the same text. It seems that they all use 
the same source - and sometimes the writer did not even copy correctly. What do you think about this: 



(from a master degree thesis!) 


Nice binomial coefficients ! !?? 
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For creating an 5-dimensional Faure sequence one has to start with a Halton sequence based on the 
smallest prime number p > s. These are the first dimension coordinates of the 5-dimensional points. 
We receive the coordinates of the next dimension by multiplying the p-based vectors - which finally 
give the decimal values - recursively by a matrix containing binomial coefficients mod p. 

It’s better to demonstrate this: 


We will generate a 4-dimensional Faure sequence. So we start with a Halton sequence based onp = 5. 
Take the 100 th element: 100 = 400 5 , which is reflected 0.004 5 giving the decimal number 0.032. 


conv_rev(100 , 5) = [0, 0, 4] 


(ha"lton_seq(100 , 
(ha"lton_5eq(100 , 



4 

5)) 

H 

100 

125 

5)) 

.d0. 032 

100 



We need the vector [0, 0, 4] as column vector which must be multiplied recursively by an upper 
triangular matrix of binomial coefficients with appropriate dimension: 


f r 0 r \ r-' i2 

C'n C'n C-'n C'r 


0 

0 

0 

0 


c c c, 


0 

0 

0 


0 

0 


C; C 


2 

Cl 

0 


mod p. This is in our case: 


■■J 


Y oA 


w 


C i\ 


/i\ 


vly 


2Y 

vOy 
f A 




\ 2 jj 


o 


2 

1 


So let’s perform the repeated multiplication mod p. DERIVE ’s ITERATES-function is an excellent 
tool for jobs like this: 




' 0 ' 

> 



0 


4 


1 


1 


ITERATES 

M0D(trm(3)« v, 5), v, 

0 

, 3 

= 


0 

f 

3 

f 

1 

P 

4 



* 

. 4 . 

J 



4 


4 


4 


4 


' 4 

4 3 4 1 


1 


4 

1 


4 


4 



, + + , + + , + + 

. 125 5 25 125 5 25 125 5 25 125 . 


4 119 34 49 


125 125 125 125 

[0.032, 0.952, 0.272, 0.392] 

Point #100 has coordinates (0.32, 0.952, 0.272, 0.392). Later we will have the occasion to check if this 
calculation is correct. 

I start implementing the algorithm in DERIVE defining two auxiliary functions: 

conv_rev(njmber , base, c) := 

Prog 
c := [] 

Loop 

#1: If number < base exit 

c := APPEND([MGD(number , base)], c) 
number := FLOORfnumber, base) 

REVER5E(APPEND([number] , c)) 

#2: trm(k) := VECTOR (VECTOR (COME (l , j), i, 0, k - 1), j, 0, k - 1) 
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con_rev(number , base) converts and reflects, giving a vector and t rm(k) generates the appropriate 
upper trianglular matrix of binomial coefficients. 

fanre(n, s, p, ir, r, fan) := 

Prog 

p := IF(PRINIE?Cs) p s, NEXT_PRIME(s )) 
fan := [] 
ir := 1 
Loop 

If ir > n exit 
b := [conv_rev(ir p p)] 1 

r := MOD(ITERATES(trm(DIMCv)).v p v P b P s - 1) P p) 
r := [VECTORCv 1 - VECTOR (p A (-k) P k P 1 P DIM(v)) p v P r)] l il 
fan := APPENDffan, r) 
i r : + 1 
fan 

The first test: 


fan re (11, 

" 0.333333 
0.666666 
0 . 111111 
0 . 444444 
0.777777 
0 . 222222 
0.555555 
0.883388 
0 . 0370570 
0.370570 
. 0.703703 

Second test: 


0.333333 
0.666666 
0 . 444444 
0.777777 
0 . 111111 
0 . 338333 
0 . 222222 
0.555555 
0.592592 
0.925925 
0.259259 


0.333333 " 

0.666666 

0.777777 

0 . 111111 

0 . 444444 

0.555555 

0.333833 

0.222222 

0.431431 

0.814814 

0.148148 . 


Tahuf; } 

TBkEE DUtWHONAL TAIWM SEQUENCES 


1 

pa 


1 


4 ! 

* 


] 

1 

0 

0 

1/3 

10 

1/3 

2 

2 

c 

0 

2/3 

2/3 

2/3 

3 

0 

i 

0 

i ft 

4/9 

7ft 

4 

1 

i 

Q. 

Aft 

1 1/9 

1ft 

5 

2 

i 

0 

7ft 

W 

4ft 

& 

0 

7 

0 

7ft 

i m 

5ft 

7 

1 

2 

0 

Sft 

2ft 

Eft 

£ 

2 

2 

0 

m 

5ft 

2ft 

9 

o 

O 

] 

1/27 

16/27 

13/17 

10 

1 

c 

] 

10/27 

23/27 

22/27 

11 

2 

0 

1 

L 9/27 

7/27 

4/27 
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fan re (10, 6) 

" 0.1428571428 
0.2857142857 
0.4285714285 
0.5714285714 
0.7142857142 
0.8571428571 
0.02040816326 
0.1632653061 
0.3061224489 
. 0.4489795918 


0.1428571428 

0.2857142857 

0.4285714285 

0.5714285714 

0.7142857142 

0.8571428571 

0.1632653061 

0.3061224489 

0.4489795918 

0.5918367346 


0.1428571428 

0.2857142857 

0.4285714285 

0.5714285714 

0.7142857142 

0.8571428571 

0.3061224489 

0.4489795918 

0.5918367346 

0.7346938775 


0.1428571428 

0.2857142857 

0.4285714285 

0.5714285714 

0.7142857142 

0.8571428571 

0.4489795918 

0.5918367346 

0.7346938775 

0.8775510204 


0.1428571428 

0.2857142857 

0.4285714285 

0.5714285714 

0.7142857142 

0.8571428571 

0.5918367346 

0.7346938775 

0.8775510204 

0.02040816326 


0.1428571428 " 

0.2857142857 

0.4285714285 

0.5714285714 

0.7142857142 

0.8571428571 

0.7346938775 

0.8775510204 

0.02040816326 

0.1632653061 . 
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For my second test I use an Excel-worksheet which can be downloaded from [11]. 
It is an LDP (= Low Discrepancy Points) Generator. 


Iterations 

100 



LDP Generation Runtime 



Dimensions 

6 



Start 

18.09.2013 10:09 




Prime 


7 



Finish 

13.09.2013 10:09 









RunTime 

00:00:04 




































Output - Low Discrepancy Points 





















0.142857 

0.142357 

0.142357 

0.142357 

0.142357 

0.142357 




0.285714 

0.235714 

0.285714 

0.235714 

0.235714 

0.235714 




0.428571 

0.423571 

0.423571 

0.423571 

0.423571 

0.423571 




0.571428 

0.571429 

0.571429 

0.571429 

0.571429 

0.571429 




0.714236 

0.714236 

0.714236 

0.714236 

0.714286 

0.714236 




0.357143 

0.357143 

0.357143 

0.357143 

0.357143 

0.357143 




0.020403 

0.163265 

0.306122 

0.44893 

0.591337 

0.734694 




0.163265 

0.306122 

0.44393 

0.591337 

0.734694 

0.377551 




0.306122 

0.44393 

0.591337 

0.734694 

0.377551 

0.020403 




0.44898 

0.591337 

0.734694 

0.377551 

0.020403 

0.163265 




0.591337 

0.734694 

0.377551 

0.020403 

0.163265 

0.306122 




0.734694 

0.377551 

0.020403 

0.163265 

0.306122 

0.44893 




0.377551 

0.020403 

0.163265 

0.306122 

0.44393 

0.591337 




0.040316 

0.326531 

0.612245 

0.397959 

0.133673 

0.469333 





Please compare! 

I use LPD for calculating the 100 th point in 4D-space from above: 


99 

0.984 

0.224 

0.864 

0.704 

100 

0.032 

0.952 

0.272 

0.392 

101 

0.232 

0.152 

0.472 

0.592 


It seems that my program works properly. 

I finish this contribution demonstrating the 
distribution of 500 points of 2D-Faure points. The 
first 125 points are the red ones. 



You can find a numerical approximation of a multiple integral on the next page. 
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QMC with TI-NspireCAS 



halton 2000, [o o], 


7T 

- 2 

2 


3,5 


Coordinates in lists halx S haly" 


Per tig 


count! f(seq [/( h a /x[ A;] ) > h a fy[ A;] J A; 1 , d im( h a lx ) ), ? 


? =tme 


1631 


1631-tt 


2.56197 


2000 


faure{ 10,6) 


QMC-numbers in columns of matrix fau 


Per tig 


fau 


0.142857 

0.142857 

0.142857 

0.142857 

0.142857 

0.14285': 

0.285714 

0.285714 

0.285714 

0.285714 

0.285714 

0.28571' 

0.428571 

0.428571 

0.428571 

0.428571 

0.428571 

0.42857! 

0.571429 

0.571429 

0.571429 

0.571429 

0.571429 

0.57142< 

0.714286 

0.714286 

0.714286 

0.714286 

0.714286 

0.71428( ( 

0.857143 

0.857143 

0.857143 

0.857143 

0.857143 

0.85714: 

0.020408 

0.163265 

0.306122 

0.44898 

0.591837 

0.73469^ 

0.163265 

0.306122 

0.44898 

0.591837 

0.734694 

0.87755! 

0.306122 

0.44898 

0.591837 

0.734694 

0.877551 

0.02040^ 

0.44898 

0.591837 

0.734694 

0.877551 

0.020408 

0.16326? 


■faure" erfolgreich gespeichert 


^))) t 


Define faure(n J s)= 

Prgm 

Local Pfir^e 
p:=nextprime[s) 

/au:=newMat(ff,s) 
ir = 1 

While ir<n 
b : = ( 1 ist> mat( c o rtv_r 
For ?jljS 

dim(fi)[l] 

fau[ir,i]= (mod(*[A;l]j)' 

k= 1 

^:=mod(/™(dim(z?)[l])- b,p) 

EndFor 
ir=ir + 1 
EndWhile 

Disp "QMC-numbers in columns of matrix fau" 
EndPrgm 


-k' 


1 3/99 


You can download the respective quasi-random.tns-file for TI-Nspire CAS. 
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Remaining Challenges 

Challenge #1: Add some introductory examples how QMC-methods are used in financial 
mathematics (Black-Scholes, option pricing, . . .). In many applications the quasi random numbers 
must be transformed to normal distributed random numbers in order to simulate a Brownian motion, 

Challenge #2: Sobol’ sequences are more complicated. Try to program the respective algorithm [9] . 

References 

There are so many websites (my “can of worms”) that I can only give a selection: 


[1] 

[2] 

http://planning.cs.uiuc.edu/node2 1 0.html, (Full book with many algorithms) 
http://marcoagd.usuarios.rdc.puc-rio.br/quasi mc.html, (including Excel-files) 

[3] 

http://www.actuaries.org/AFIR/Colloquia/Caims/Bovle Tan.pdf 

[4] 

http://www.math.uwaterloo.ca/~clemieux/papers/iccs.pdf (QMC in Finance) 

[5] 

http://www.smartquant.com/references/MonteCarlo/mc6.pdf (Applications of MC/QMC) 

[6] 

https://quanto.inria.fr/pdf html/mc quasi doc.pdf 

[7] 

http://www.mcqmc2012.unsw.edu.au/slides/MCOMC2012 Dick Tutorial.pdf 

[8] 

http://people.sc.fsu.edu/~iburkardt/rn src/faure/faure.html (MATLAB -Library) 

[9] 

http://web.maths.unsw.edu.au/~fkuo/sobol/ioe-kuo-notes.pdf 

[10] 

http://www.wpi.edu/Pubs/ETD/Available/etd-0 1131 04- 1 40925/unrestrictcd/krykova.pdf 

[11] 

http://www.casact.org/pubs/fomm/99spfomm/99spf337.pdf 


http://www.casact.org/pubs/fomm/99spfomm/ldpmaker.xls 


German websites 

[12] http://gcsc.uni-frankfurt.de/computational-finance/students/theses/Benny%20Xiang%20Li.pdf 

[13] http://kainhofer.com/Papers/Kainhofer Talk QMC Wissenswertes 20-06-2005.pdf 

[14] http://itp.tugraz.at/MML/MonteCarlo/MCIntro.pdf (Grundlagen der MC-Methode) 

[14] http://www.schulmediothek.de/fileadmin/pdf/Quasi-Monte Carlo Methode.pdf 

[15] http://pauh.uni-muenster.de/~lemm/seminarSS08/horstmann vortragqmc.pdf 
and all the respective Wikipedia-sites. 


Approximation of the area of the maple leaf applying a QMC-method: 

1024 Hammersley QMC-points give an approximation for the Scientific American maple leaf of 
41.41% of the square. The real area is 41.85% of the square 
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A Tribute to DERIVE - A Challenge for TI-NspireCAS 

When I was busy with the two ACA-contributions of our Canadian friends I got the idea that it might 
be useful trying to keep the many Utility files of DERIVE alive. They are true gems for the 
DERIVIANs and why not making them accessible for Nspire Users in form of respective libraries. To 
begin with I started with LinearAlgebra.mth. There is a library for linear algebra linalgcas.tns prod- 
uced by PhilippeFortin which is very useful but it does not contain all tools conained in the DERIVE 
file. Examples will be given in the next DNL. 

I wrote Michel Beaudin and he immediately answered proposing that we could prepare a common talk 
at TIME 2014. 1 would like to invite you all to contribute for our BIG PLAN - very well considering 
that LinearAlgebra.mth is an easy task compared with e.g. ZetaFunctions.mth. 

It would be great to hear from you! 

Row Reduction Primitives 

scale_element(v,i,s) returns a copy of vector v with the i th element multiplied by expression s. 

If v is a matrix then the i th row is multiplied by s. 

swap_elements(v,ij) returns a copy of vector v (or list v) with the i th and j th element interchanged. 

If v is a matrix then the respective rows are interchanged. 

swap_cols(A,i,j) returns a copy of matrix A (or list v) with the i th and j th column interchanged. 

This function is not part of LinearAlgebra.mth. 

subtr_elements(A,i,j\s) returns a copy of matrix A in which the i th row is replaced with the i th row 
minus s times the j th row. Is identical with Nspire's mRowAdd(-s,A,j,i). 

forceO(A,i,j\p) returns a copy of matrix A in which element Ai,j is forced to zero by substracting 
an appropriate multiple of pivot row p from row i. 

pivot(A,i,j) returns a copy of matrix A in which the elements in the j th column below the i th row are 
forced to zero 

Inverses and Adjoints by Cofactor Expansion 

delete(v,i) returns a copy of vector or list v with all elements omitted whose numbers are given in the 
2 nd parameter. If v is a matrix then the respective rows are omitted. 

minor(A,i,j) returns a copy of matrix A in which the i th row and the j th column are omitted. 

cofactor(A,i,j) returns the numerator of the row j col i element of the inverse of square matrix A. The 
corresponding element of the inverse is this numerator divided by the determinant of A. 

adjoint(A) returns adjugate of square matrix A. It is the transpose of the matrix of the cofactors of A. 

Null Space (Kernel), Eigenvalues and Eigenvectors 

The respective DERIVE-^rogram contains some auxiliary functions which are not available in Nspire 
(e.g. POSITION, FIRST, REST a.o.). They must be predefined: 

position(e,v) gives the position of first appearance of element e in vector v 

first(v) returns the first element of a vector/list or matrix (first row). 

rest(v) returns all but the first element of a vector/list or matrix (first row). 

rank(A) returns the the rank of matrix A. (from the Nspire library linalgcas.tns) 

null_space(A): The null space is the set of vectors x satisfying the matrix equation Ax=0. 

Ph. Fortin's function kernelvectors gives the answer as a row vector in form of a linear combination: 

eigenvalues, exact_evc and approx_evc provide a better functionality. 


